/**
* @file decompose.cpp
* March, 2014
* 
* @brief class to decompose quad-mesh in mapping mesh.
* 
* @author Antonio Miranda (amiranda@tecgraf.puc-rio.br )
* 
* @see 
*/

// includes
//----------

#include "decompose.h"
#include "topology.hpp"


//
//////////////////////////////////////////////////////////////////////////
decomposeQuadMesh::decomposeQuadMesh( cTopology *top )
{
  m_topol = top;
}

//
//////////////////////////////////////////////////////////////////////////
decomposeQuadMesh::~decomposeQuadMesh(void)
{
}

//
//////////////////////////////////////////////////////////////////////////
bool decomposeQuadMesh::process()
{
  int i, num_adj_edges;

  // build loops
  
  // mark all edges as 0
  int num_edges = m_topol->NumEdges();
  for (i = 0; i < num_edges; ++i)
    m_topol->SetEdgeInfo (i, 0);

  // find nodes that have adjacent number of edges different of four
  // search only in internal nodes
  int num_nodes = m_topol->NumNodes();
  vector <int> start_nodes;
  for (i = 0; i < num_nodes; ++i)
  {
    m_topol->SetNodeInfo (i, 0);

    if (m_topol->IsBoundaryNode (i)) // internal nodes is 
      continue;

    num_adj_edges = m_topol->NumAdjEdgeToNode (i);
    if (num_adj_edges != 4)
      start_nodes.push_back(i);
  }


  // initial temporary dividing curve that decompose mesh
  vector < vector<int> > temp_pts;
  for (size_t k = 0; k < start_nodes.size(); ++k)
  {

    // adjacent edges to start nodes
    num_adj_edges = m_topol->NumAdjEdgeToNode (start_nodes[k]);
    for (int m = 0; m < num_adj_edges; ++m)
    {
      vector <int> points;
      int  adj_edge = m_topol->AdjEdgeToNode (start_nodes[k], m);

      // get internal points
      getInternalIdPoints(start_nodes[k], adj_edge, points);

      if (points.size() > 0)
        temp_pts.push_back(points);
    }
  }

  // split
   vector < vector<int> > split_curve;
  for (size_t k = 0; k < temp_pts.size(); ++k)
    splitCurve (temp_pts[k], split_curve);

  // get points
  getDividingPoints(split_curve);

  getQuadConnect (split_curve);


  writeDecompToVTK("edges.vtk", split_curve);

  writeMeshToVTK ("quad_elements.vtk");


  return true;
}



//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::getInternalIdPoints( int init_pts, int init_edge, vector <int> &points )
{

  int ei, ej, info = 0;

  // insert initial point
  points.push_back(init_pts);

  // run through edges, add other points
  int curr_edge = init_edge;
  int last_pts  = init_pts;

  do
  {
    // set flag of current edges as 1
    m_topol->SetEdgeInfo(curr_edge, 1);

    // get adjacent nodes to current edge
    m_topol->AdjNodeToEdge (curr_edge, &ei, &ej);

    // add next point
    if (last_pts != ei)
    {
      points.push_back(ei);
      last_pts = ei;
    }
    else
    {
      points.push_back(ej);
      last_pts = ej;
    }

    // get next edge on curve
    curr_edge = nextOppositeEdge (curr_edge, last_pts);

    info = 0;
    if (curr_edge != -1)
      info = m_topol->GetEdgeInfo(curr_edge);

  }
  while (curr_edge != -1 && info == 0);

}


//
//////////////////////////////////////////////////////////////////////////
int decomposeQuadMesh::nextOppositeEdge(int curr_edge, int v_base )
{
  int i, elm0, elm1, elm_0, elm_1;

  // get adjacent elements to edge
  if (m_topol->NumAdjElemToEdge  (curr_edge) != 2)
    return -1;
  elm0 = m_topol->AdjElemToEdge (curr_edge, 0);
  elm1 = m_topol->AdjElemToEdge (curr_edge, 1);


  int n = m_topol->NumAdjEdgeToNode (v_base);
  for (i = 0; i < n; ++i)
  {
    int curr_edge = m_topol->AdjEdgeToNode (v_base, i);

    if (m_topol->NumAdjElemToEdge (curr_edge) != 2)
      return -1;

    elm_0 = m_topol->AdjElemToEdge (curr_edge, 0);
    elm_1 = m_topol->AdjElemToEdge (curr_edge, 1);

    if (elm_0 != elm0 && elm_0 != elm1 && 
      elm_1 != elm0 && elm_1 != elm1 ) 
      return curr_edge;
  }

  return -1;
}




//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::
splitCurve( vector <int> &temp_curve, vector < vector<int> > &div_curve )
{
  int curr_node, n_adj_edges;
  vector <int> new_curve;


  // init a new curve with the initial point
  new_curve.push_back(temp_curve[0]);

  // for all points in dividing curve
  size_t n_pts = temp_curve.size();
  for (size_t k = 1; k < n_pts; ++k)
  {
    // current node
    curr_node = temp_curve[k];
    new_curve.push_back(curr_node);

    // number of adjacent edges to node
    n_adj_edges = m_topol->NumAdjEdgeToNode (curr_node);

    // check if all adjacent edges to node is marked, so this point divide 
    // the curve
    bool div_point = true;
    for (int m = 0; m < n_adj_edges; ++m)
    {
      int  adj_edge = m_topol->AdjEdgeToNode (curr_node, m);

      // get internal points
      if (m_topol->GetEdgeInfo(adj_edge) == 0)
      {
        div_point = false;
        break;
      }
    }

    // end a dividing curve and star another
    if (div_point)
    {
      div_curve.push_back (new_curve);
      
      new_curve.clear();
      new_curve.push_back(curr_node);
    }
  }
  
}

//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::getDividingPoints( vector < vector<int> > &div_curve )
{
  int index = 0;
  double coord[3];

  // set all info-flag node = -1;
  for (size_t k = 0; k < div_curve.size(); ++k)
  {
    // get first and last point in a curve
    int id_test[2] = {div_curve[k][0], div_curve[k][div_curve[k].size()-1]};
    
    for (int m = 0; m < 2; ++m)
      m_topol->SetNodeInfo (id_test[m], -1);
  }

  // set info-flag with new ids
  int idx_node = 0;
  for (size_t k = 0; k < div_curve.size(); ++k)
  {
    int id_test[2] = {div_curve[k][0], div_curve[k][div_curve[k].size()-1]};

    // test in the point exist
    for (int m = 0; m < 2; ++m)
    {
      int info = m_topol->GetNodeInfo (id_test[m]);

      if (info == -1)
      {
        m_topol->SetNodeInfo (id_test[m], idx_node);
        m_topol->GetCoordNode(id_test[m], coord);
        m_quad_pts.push_back(coord[0]);
        m_quad_pts.push_back(coord[1]);
        m_quad_pts.push_back(coord[2]);
        ++idx_node;
      }
    }
  }

}


//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::getQuadConnect (vector < vector<int> > &div_curve)
{

  // set all adjacent edges adjacent to initial and final node at
  // dividing curve with info-flag = id_curve + 1
  for (size_t k = 0; k < div_curve.size(); ++k)
  {
    int edge_ids[2][2] = { {div_curve[k][0], div_curve[k][1]},
                          {div_curve[k][div_curve[k].size()-2], div_curve[k][div_curve[k].size()-1]} } ;
    // 
    for (int m = 0; m < 2; ++m)
    {
      int curr_edge = m_topol->GetEdge (edge_ids[m][0], edge_ids[m][1]);
      m_topol->SetEdgeInfo (curr_edge, k+1);
    }
  }


  // set all adjacent elements no starting node of curve as info-flag = -1
  for (size_t k = 0; k < div_curve.size(); ++k)
  {
    int ne, *elem;
    m_topol->AdjElemToNode (div_curve[k][0], &ne, &elem);
    for (int m = 0; m < ne; ++m)
      m_topol->SetElemInfo(elem[m], -1);
    delete []elem;
  }

  // obtain elements starting from curves
  int idx_elem = 0;
  for (size_t k = 0; k < div_curve.size(); ++k)
  {
    // adjacent element to edge in sequence j -> j
    int edge_i = div_curve[k][0];
    int edge_j = div_curve[k][1];

    int adj_elem = adjElemToEdgeIJ (edge_j, edge_i);
    
    // do not loop if the element was visited
    if (m_topol->GetElemInfo(adj_elem) != -1)
      continue;

    int start_node, curr_node = -1;
    start_node = div_curve[k][0];

    // init new element
    m_quad_conn.push_back(4);
    m_quad_conn.push_back(m_topol->GetNodeInfo(start_node));

    while (true)
    {
      // id edge with nodes i and j
      int curr_edge  = m_topol->GetEdge (edge_i, edge_j);
      int curr_curve = m_topol->GetEdgeInfo (curr_edge) - 1;  // must decrease 1 because we add +1 in beginning.

      // end edge in current curve
      if (edge_i == div_curve[curr_curve][0])
      {
        edge_i = div_curve[curr_curve][div_curve[curr_curve].size()-2];
        edge_j = div_curve[curr_curve][div_curve[curr_curve].size()-1];
      }
      else
      {
        edge_i = div_curve[curr_curve][1];
        edge_j = div_curve[curr_curve][0];
      }


      // check of final loop
      if (start_node == edge_j)
        break;
      else
        m_quad_conn.push_back(m_topol->GetNodeInfo(edge_j));

      // adjacent element to edge in sequence j -> i, mark with zero
      adj_elem = adjElemToEdgeIJ (edge_j, edge_i);
      m_topol->SetElemInfo(adj_elem, 0);
      

      // next edge in other curve
      edge_i = edge_j;
      edge_j = prevNodeElem (adj_elem, edge_i);
    }
  }
  
}



//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::getDMesh( int *np, double **pts, int *ne, int **conn )
{

}

//
//////////////////////////////////////////////////////////////////////////
cTopology * decomposeQuadMesh::getDMesh()
{
  cTopology *quad_mesh = NULL;

  return quad_mesh;
}

//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::
nextEdgeOnQuad( int edge_i, int edge_j, int nxt_edge_i, int nxt_edge_j )
{
  // current edge
  int curr_edge = m_topol->GetEdge(edge_i, edge_j);
  
  // get id dividing curve
  int id_curve = m_topol->GetEdgeInfo(curr_edge);

}

//
//////////////////////////////////////////////////////////////////////////
int decomposeQuadMesh::adjElemToEdgeIJ( int edge_i, int edge_j )
{
  int k, m, nxt;
  // current edge
  int curr_edge = m_topol->GetEdge(edge_i, edge_j);

  // adjacent elements to edge
  int n_adj_elem = m_topol->NumAdjElemToEdge (curr_edge);
  for (k = 0; k < n_adj_elem; ++k)
  {
    int curr_elem = m_topol->AdjElemToEdge(curr_edge, k);
    int n_nodes   = m_topol->GetElemNNodes (curr_elem);
    int *conn     = m_topol->GetElemConn(curr_elem);
    
    for (m = 0; m < n_nodes; ++m)
    {
      nxt = (m+1)%n_nodes;
      if (conn[m] == edge_i && conn[nxt] == edge_j)
        return curr_elem;
    }
  }

  return -1;
}

//
//////////////////////////////////////////////////////////////////////////
int decomposeQuadMesh::prevNodeElem( int curr_elem, int id_node )
{
  int m, nxt;

  int n_nodes   = m_topol->GetElemNNodes (curr_elem);
  int *conn     = m_topol->GetElemConn(curr_elem);

  for (m = 0; m < n_nodes; ++m)
  {
    nxt = (m+1)%n_nodes;
    if (conn[nxt] == id_node)
      return conn[m];
  }

  return -1;
}

//
//////////////////////////////////////////////////////////////////////////
void decomposeQuadMesh::
writeDecompToVTK( string filename, vector < vector<int> > ids_pts )
{
  FILE *arquivo = NULL;

  arquivo = fopen(filename.c_str(),"w");
  if (arquivo == NULL)
    return;

  // number of points
  size_t npts = 0, nedges = 0;
  for (size_t k = 0; k < ids_pts.size(); ++k)
  {
    npts   += ids_pts[k].size();
    nedges += (ids_pts[k].size() - 1);
  }

  fprintf(arquivo,"# vtk DataFile Version 2.0\n");
  fprintf(arquivo,"Decompose unstructure quad-mesh in mapping mesh\n");
  fprintf(arquivo,"ASCII\n\n");

  fprintf(arquivo,"DATASET UNSTRUCTURED_GRID\n");
  fprintf(arquivo,"POINTS %d float\n", (int) npts);

  //////////////////////////////////////////////////////////////////////////
  double coord[3];
  for (size_t k = 0; k < ids_pts.size(); ++k)
  {
    for (size_t m = 0; m < ids_pts[k].size(); ++m)
    {
      m_topol->GetCoordNode(ids_pts[k][m], coord);
      fprintf(arquivo,"%f    %f    %f\n", coord[0], coord[1], coord[2]);
    }
  }


  int base_id_num = 0;

  fprintf(arquivo,"CELLS %d %d\n", nedges + (int)(ids_pts.size()*2), nedges*3 + (int)(ids_pts.size()*4));
  for (size_t k = 0; k < ids_pts.size(); ++k)
  {

    // write initial vertex
    fprintf(arquivo,"1    %d\n", base_id_num);

    for (size_t m = 0; m < ids_pts[k].size()-1; ++m)
      fprintf(arquivo,"2    %d    %d\n", base_id_num+m, base_id_num+m+1);

    base_id_num += (int) ids_pts[k].size();

    // write end vertex
    fprintf(arquivo,"1    %d\n", base_id_num-1);
  }


  fprintf(arquivo,"CELL_TYPES %d\n", nedges + (int)(ids_pts.size()*2));
  for (size_t k = 0; k < ids_pts.size(); ++k)
  {
    // write vertex
    fprintf(arquivo,"1\n");

    for (size_t m = 0; m < ids_pts[k].size()-1; ++m)
      fprintf(arquivo,"3\n");

    // write vertex
    fprintf(arquivo,"1\n");
  }


  fclose(arquivo);
}

void decomposeQuadMesh::writeMeshToVTK( string filename )
{
  FILE *arquivo = NULL;

  arquivo = fopen(filename.c_str(),"w");
  if (arquivo == NULL)
    return;

  // number of points
  int npts = (int) m_quad_pts.size() / 3;

  fprintf(arquivo,"# vtk DataFile Version 2.0\n");
  fprintf(arquivo,"Decompose unstructure quad-mesh in mapping mesh\n");
  fprintf(arquivo,"ASCII\n\n");

  fprintf(arquivo,"DATASET UNSTRUCTURED_GRID\n");
  fprintf(arquivo,"POINTS %d float\n", (int) npts);

  //////////////////////////////////////////////////////////////////////////
  for(int i=0; i<npts; i++)
    fprintf(arquivo,"%f    %f    %f\n", m_quad_pts[i*3+0], m_quad_pts[i*3+1], m_quad_pts[i*3+2]);


  int NumElem = (int) m_quad_conn.size() / 5;
  fprintf(arquivo,"CELLS %d %d\n", NumElem, NumElem*5);
  for(int i=0; i<NumElem; i++)
  {
    fprintf(arquivo,"%d", 4);
    for(int j=0; j< 4; j++)
      fprintf(arquivo,"  %d", m_quad_conn[i*5+j+1]);
    fprintf(arquivo,"\n");
  }

  fprintf(arquivo,"CELL_TYPES %d\n", NumElem);
  for(int i=0; i<NumElem; i++)
      fprintf(arquivo,"9\n");

  fclose(arquivo);
}




